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When atoms and molecules are irradiated by an x-ray free-electron laser (XFEL), they are highly 
ionized via a sequence of one-photon ionization and relaxation processes. To describe the ionization 
dynamics during XFEL pulses, a rate equation model has been employed. Even though this model 
is straightforward for the case of light atoms, it generates a huge number of coupled rate equations 
for heavy atoms like xenon, which are not trivial to solve directly. Here, we employ the Monte Carlo 
method to address this problem and we investigate ionization dynamics of xenon atoms induced 
by XFEL pulses at a photon energy of 4500 eV. Charge state distributions, photo-/ Auger electron 
spectra, and fluorescence spectra are presented for x-ray fluences of up to 10 13 photons//xm 2 . With 
the photon energy of 4500 eV, xenon atoms can be ionized up to +44 through multiphoton absorption 
characterized by sequential one-photon single-electron interactions. 

PACS numbers: 32.80.Fb, 32.90.+a, 41.60.Cr, 02.70.Uu 



I. INTRODUCTION 

The recent advent of x-ray free-electron lasers 
(XFEL) [1-3] enables us to explore new frontiers of sci- 
ence [4], for example, femtosecond x-ray imaging [5-9] 
and warm dense matter [10]. A series of experiments con- 
ducted at the Linac Coherent Light Source (LCLS) [11] 
have shown how ultraintense and ultrashort x rays in- 
teract with various systems: light atom (Ne) [12, 13], 
molecule (N2) [14-16], heavy atom (Xe) [17], and solid 
(Al) [18]. ^ 

The ionizing XFEL-matter interaction is one of the 
most fundamental processes that affects all XFEL appli- 
cations. As demonstrated theoretically [19] and experi- 
mentally [12], the electronic response to an XFEL pulse 
is characterized by a sequence of one-photon ionization 
and relaxation events. In the x-ray regime, photoab- 
sorption predominantly ionizes an inner-shell electron. 
The resulting inner-shell vacancy is filled via radiative 
(fluorescence) and/or non-radiative (Auger and Coster- 
Kronig) transitions. Then the extremely large number 
of x-ray photons within an ultrashort XFEL pulse can 
keep ionizing after or even before these relaxation pro- 
cesses are over [19, 20]. As a result, atoms or molecules 
become highly ionized after absorbing several photons 
sequentially. To describe ionization dynamics, we em- 
ploy a rate equation model, which demonstrates good 
agreement with experiments conducted at LCLS [12, 13]. 
Tracking populations via rate equations is sufficient to 
describe the ionization dynamics during XFEL pulses, 
mainly because the coherence time of current XFEL 
sources is much shorter than the time scale of popu- 
lation changes. The ionization dynamics in an XFEL 
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pulse differ from those at a third-generation x-ray syn- 
chrotron radiation source, where one-photon absorption 
is dominant, and from multiphoton strong-field ioniza- 
tion, where many photons are simultaneously absorbed 
to ionize a single electron. Understanding radiation dam- 
age mechanisms [21] including ionization dynamics is of 
central importance for single-shot imaging of individual 
molecules [22-24]. 

To probe ionization dynamics induced by an XFEL, 
one can collect all particles generated in the interac- 
tion between XFEL and matter. Photoionization and 
Auger (Coster-Kronig) decay produce electrons, and flu- 
orescence produces photons. Also highly charged ions 
are generated via multiphoton multiple ionization. It is 
possible to simultaneously measure all those particles by 
means of the CFEL-ASG Multi-Purpose (CAMP) instru- 
ment [25] , which has been successfully applied to a study 
of the XFEL-heavy atom interaction [17] and to single- 
shot imaging experiments [26-30]. 

Ionization of heavy atoms irradiated by XFEL pulses 
has attracted considerable attention, not only because 
the heavy atom has many electrons to be ionized but 
also because it has a rich manifold of ionization channels 
involving complex inner-shell decay cascades. A recent 
study proposes a resonance-enabled x-ray multiple ion- 
ization mechanism for heavy atoms to reach high charge 
states beyond those expected from the straightforward 
sequential ionization model [17]. Also, ionization dynam- 
ics of heavy atoms embedded in macromolecules deliver 
a novel way to determine macromolecular structure from 
femtosecond nanocrystallography data using XFEL [31]. 
The theoretical treatment of the XFEL-heavy atom in- 
teraction is challenging because its dynamics involve a 
huge number of possible pathways. For example, the 
simple rate equation model for Xe M-shell ionization re- 
quires more than one million coupled rate equations [17]. 
This number becomes even much larger when resonantly 
excited Rydberg states are taken into account. There- 
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fore, it is necessary to develop a computational tool that 
can handle heavy atoms and provide full information on 
ions, electrons, and photons. 

In the present work, we employ the Monte Carlo 
method to solve the large number of coupled rate equa- 
tions. There have been extensive studies using the Monte 
Carlo procedure for decay pathways of an inner-shell va- 
cancy produced by x-ray synchrotron or electron cap- 
ture [32-44]. In particular, the inner-shell decay pro- 
cess for heavy atoms such as iodine and xenon [34, 36- 
38] has brought much interest because of their relevance 
for medical applications [38, 45]. In conventional imple- 
mentations of the Monte Carlo method, atomic data for 
multiple-hole ions are usually assumed to be the same as 
for the singly-ionized atom [36] or they are scaled from 
the singly-ionized atom according to the number of va- 
lence electrons [46]. To the best of our knowledge, no 
Monte Carlo calculation has been done using atomic data 
that are individually calculated for all possible multiple- 
hole configurations. In the current implementation, we 
perform electronic structure calculations individually for 
all possible configurations in order to obtain a whole set 
of atomic data. This is important for ionization dynam- 
ics in XFEL radiation because of the production of a 
broad range of charge states. The higher the ionic charge, 
the greater are the deviations from the singly-ionized 
atom. As the charge state increases, some Auger tran- 
sitions become energetically forbidden [39]. Thus, the 
detailed electronic structure for each configuration mat- 
ters to atomic data calculations and eventually to ion- 
ization dynamics simulations. Another distinct aspect of 
the Monte Carlo implementation described in this paper 
is the availability of temporal information on electronic 
dynamics. In ordinary Monte Carlo simulations of elec- 
tronic decay cascades, the time variable has not been of 
main concern, and in many cases decay rates are simply 
given by relative rates. However, it is important to know 
details of ionization dynamics, especially when their time 
scale is comparable with the time scale of nuclear dynam- 
ics [47]. 

The structure of this paper is as follows. In Sec. II, we 
describe the rate equation model for ionization dynamics 
and present a Monte Carlo implementation for solving a 
large set of coupled rate equations. In Sec. Ill, we present 
atomic data of Xe and time-dependent ionization path- 
ways from Monte Carlo simulations. We discuss ioniza- 
tion dynamics of Xe in intense x-ray pulses at a photon 
energy of 4500 eV by analyzing ion, electron, and photon 
spectra. We conclude with a summary and outlook in 
Sec. IV. 



II. THEORY AND NUMERICAL DETAILS 

A. Ionization dynamics 

To simulate ionization dynamics in intense x-ray 
pulses, we employ a rate equation approach based on se- 



quential one-photon ionization and relaxation steps. This 
rate equation model was introduced in connection with 
XFEL-atom interactions by Rohringer and Santra [19], 
extended to x-ray scattering dynamics [24], generalized 
to arbitrary elements [48], and has been successfully ap- 
plied to explain recent LCLS experiments [12, 13, 17]. 
A rate equation model for non-local thermal equilibrium 
plasma [49, 50] has been applied to a warm dense matter 
study at LCLS [18]. 

Here, we summarize procedures underlying the rate 
equation model. For interested readers, theoretical back- 
ground [51] and detailed descriptions [24] are available. 
For a given atom, we construct all possible electronic con- 
figurations {/} that may be formed by removing zero, 
one, or more electrons, from the neutral ground config- 
uration. The orbital structures are optimized with the 
Hartree-Fock-Slater (HFS) method for each configura- 
tion. We include all possible one-photon ionization and 
relaxation processes for each configuration, i.e., the sub- 
shell photoionization cross sections for a given photon 
energy, Auger (Coster-Kronig) decay rates, and fluores- 
cence rates are calculated for every single configuration. 
We calculate shake-off branching ratios, based on the 
sudden change approximation [52]. In the present work, 
shake-off processes are included for all photo-induced 
processes of neutral Xe. The calculated cross sections 
and rates serve as input parameters for a set of rate equa- 
tions of the form, 

, all config. 

jm)= £ [iwp^-iwp^)], (i) 

where Pi is the population of the 7th configuration, and 
r/_j./< is the rate for a transition from the configuration 
I to the configuration 

For the heavy atom case, the numbers of configura- 
tions and processes involved in ionization dynamics are 
very large. For example, an x-ray photon of 4500 eV 
can ionize the M-, N-, and O-shells of Xe. The num- 
ber of all possible electronic configurations constructed 
from these ionizations is 1,120,581, which is equal to 
the number of coupled rate equations to be solved, and 
the number of all possible processes under consideration 
is 43,221,650. Therefore we need to propagate in time 
a matrix of approximately one million by one million 
with ~40 million nonzero elements, in order to simulate, 
within the rate equation model, ionization dynamics of 
Xe exposed to 4500-cV XFEL pulses. To avoid a direct 
time-propagation solution of this huge matrix, we develop 
a Monte Carlo approach that efficiently solves the rate 
equations, as an extension of the xatom toolkit [48] . 



B. Monte Carlo implementation 

Our Monte Carlo description of ionization dynamics 
may be summarized as follows. For one realization (a 
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Monte Carlo trajectory), a given atom undergoes a se- 
quence of photoionization and relaxation events during 
the time propagation and eventually ends up in a final 
charge state. There are many pathways to reach the same 
final charge state. Each pathway consists of many steps 
of photoionization and relaxation, and those steps are 
stochastically determined. Probabilities are calculated 
at a given time and the time increment is optimally de- 
termined during the time propagation. If it is no longer 
possible to proceed to a further process, the time prop- 
agation ends and this trajectory is complete. We repeat 
this procedure for many trajectories to form a sufficiently 
large statistical ensemble. 

In the direct time-propagation solution, configuration 
populations are given by a fractional number represent- 
ing the probabilities for all individual configurations. On 
the other hand, for each realization of the Monte Carlo 
implementation, configuration populations are given by 
either zero or one, thus following a specific pathway of 
configuration changes induced by photoionization and re- 
laxation events. After running many trajectories, we ob- 
tain ensemble-averaged configuration populations, which 
are ideally the same as the configuration populations ob- 
tained by the direct solution. 

Here is a more detailed description of our Monte Carlo 
implementation. 

(a) Choose an initial value of At. This is also used for 
the maximum value of At. 

(b) Set up an initial configuration /, which is usually 
given by the ground configuration of a neutral atom. 
Set an initial value for the time t. 

(c) Calculate transition probabilities, {pk} for 1 < fc < 
Ap roc - Here, fc indicates an index for the transition 
process of I —> I', and N ploc is the number of all 
possible processes from /. 

!(T^J(t)At for photoionization, 
r A At for Auger (Coster-Kronig) decay, 
r F At for fluorescence, 

where J(t) is the photon flux of the x-ray pulse at a 
given time t, <t p is the photoionization cross section, 
r A is the Auger (Coster-Kronig) rate, and T F is the 
fluorescence rate. 

(d) Construct a table of processes, {Tk}. T = and 

T k = Sfc'=iPfc' for 1 < fc < A proc . Here, T Nproc gives 
the total probability to proceed to one of the N ploc 
processes, whereas 1 — Tn pioc gives the probability to 
remain in configuration /. 

(e) Choose a random number, r € [0,1]. If Tk-\ < 
r < Tk, then it proceeds to the fcth process and the 
new configuration becomes /'. For Auger (Coster- 
Kronig) decay and photoionization processes, the 
electron count corresponding to its kinetic energy bin 



is increased by one. For fluorescence, the photon 
count corresponding to its emitted photon energy bin 
is increased by one. 

(f) Adjust At according to the calculated 2jv proc . It must 
satisfy T/v -C 1 and must not be larger than the 
initial At. Then increase the time variable t by At. 

(g) Repeat (c)-(f) as long as T Nproc £ 0. If T Nproc = 0, 
there will be no further process. This finishes one 
Monte Carlo trajectory of the time-propagation cal- 
culation. The final charge state count is increased by 
one. 

(h) Run many trajectories until the results are converged. 

(i) The counts of charge state, electron energy, and pho- 
ton energy are divided by the number of trajectories. 
These histograms correspond to the ion, electron, and 
photon spectra, respectively. 

We use the following computational parameters. For 
(a), the initial and maximum value of At is 10 attosec- 
onds. For (f), At is chosen such that TW c =0.1 during 
the time propagation. For the convergence criterion of 
(h), we check charge-state populations for every 100 tra- 
jectories. In practice, 10,000 to 30,000 trajectories are 
carried out to obtain convergency of 10~ 4 for all charge- 
state populations. These Monte Carlo results fully agree 
with the direct solution to within an accuracy of 10~ 3 dis- 
crepancies. There is a tremendous reduction in the com- 
putational time. The direct solution takes about 1,760 
minutes for a Xe-atom ionization dynamics calculation 
with 16,000 time steps on the lab workstation, whereas 
the Monte Carlo implementation takes only 7 minutes 
with about 20,000 trajectories to get converged results. 



III. RESULTS AND DISCUSSION 

A. Atomic data 

Figure 1 shows orbital binding energies of the ground 
configuration of Xe 9+ as a function of the charge state 
+q. The red lines (M-shell) from the bottom correspond 
to the 3s, 3p, and 3d subshells. The green lines (iV-shell) 
correspond to the 4s, Ap, and Ad subshells, and the blue 
lines (O-shell) correspond to the 5s and hp subshells. The 
dots with circles, triangles, and rectangles, indicate that 
the corresponding subshell is fully or partially filled with 
electrons for given charge states. Here are some examples 
of the ground configuration of Xe 9+ , 

Xe 0+ : ls 2 2s 2 2p 6 3s 2 3p 6 M w As 2 Ap 6 Ad 10 5s 2 5p 6 , 
Xe 8+ : ls 2 2s 2 2p 6 3s 2 3p 6 3d 10 As 2 Ap 6 Ad 10 , 

Xe 26+ : ls 2 2s 2 2/3s 2 3/3d 10 , 

Xe 44+ : ls 2 2s 2 2/. 
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FIG. 1. (Color online) Orbital binding energies of the ground 
configuration of Xe and its charge states. The symbols (circle, 
triangle, and square) indicate that the corresponding subshell 
contains at least one electron. 



As shown in Fig. 1, the photon energy of 4500 eV is well 
above all ionization potentials of M-, N-, and O-shell 
electrons for all charge states of Xe. No resonance tran- 
sition is expected with this photon energy. Thus one can 
expect Xe 44+ as the maximum charge state after multi- 
photon multiple ionization by x rays of 4500 eV, if the 
x-ray photon fluence is high enough to remove all n>3 
electrons via a sequence of photoionization and relaxation 
processes. 

Table I compares fluorescence, Auger, and Coster- 
Kronig rates with semi-empirical calculations [53] for M- 
shell single-hole configurations of Xe. Semi-empirical cal- 
culations employ transition energies from experiments, 
whereas the present method computes them from HFS 
orbital energies. For this reason, most semi-empirical 
calculations consider only single-hole or double-hole con- 
figurations. However, we emphasize that for the present 
Monte Carlo simulations all multiple-hole configurations 
are individually calculated with the HFS method. Even 
though the rates presented in this table are summed over 
subshells X and Y, all transitions to individual subshells 
are calculated for the present Monte Carlo simulations 
and those numbers are comparable with the extensive 
table [54] for the single-hole configurations. The semi- 
empirical calculations include relativistic changes in the 
transition energies via jj-coupling [55]. On the other 
hand, the present calculations based on the nonrelativis- 
tic HFS method do not include fine-structure splittings 
(for example, between M2 and M3, or between M4 and 
M5). Without relativity, M2-M3X Coster-Kronig tran- 
sitions are energetically impossible [56, 57], so this tran- 
sition is completely absent in the present method as 
shown in Table I. In spite of these limitations of the HFS 
method, this comparison shows good agreement. 

The photoionization cross sections for neutral Xe 
are compared between the relativistic method and the 
present method in Table II. The relativistic results are 
based on the Dirac-Fock-Slater method [58] and summed 
over different total angular momentum j to make a com- 



TABLE I. Comparison of fluorescence (Mi-X), Auger (Mi- 
XY), and Coster-Kronig (Mi-MjX) rates for M-shell single- 
hole configurations of Xe. The rates are given in atomic units. 
SE represents semi-empirical calculations [53]. M23 and M45 
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TABLE II. Comparison of photoionization cross sections (in 
kilobarns) for neutral Xe. The Dirac-Fock-Slater (DFS) re- 
sults [58] are calculated at 4509 eV and the present results 
are at 4500 eV. The DFS results are summed over different 
total angular momentum j. 



Subshell DFS [58] Present 

3s 7T99 YM 

3p 25.33 24.20 

3d 16.37 16.22 

4s 1.84 1.82 

4p 4.96 4.73 

4d 2.58 2.57 

5s 0.28 0.27 

5p 0.54 0.52 



parison with the nonrelativistic case. For this case, the 
present results are in excellent agreement with the rela- 
tivistic results. 

Finally, Table III lists decay widths of single-hole con- 
figurations of Xe. The widths are calculated by the sum 
of fluorescence, Auger, and Coster-Kronig rates. The 
present results are compared with semi-empirical calcu- 
lations [53, 59, 60] and recommended values from various 
experiments and theories [61]. Bearing in mind some lim- 
itations of the nonrelativistic treatment for heavy atoms, 
the present results are in fair agreement with other avail- 
able values. 



B. Ionization pathways 

By analyzing Monte Carlo trajectories, one can re- 
trieve useful information on ionization pathways dur- 
ing XFEL pulses. Figure 2 shows 100 exemplary tra- 
jectories that are randomly chosen out of 22,200 tra- 
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TABLE III. Comparison of the decay widths (in eV) of single- 
hole configurations of Xe. EXP represents recommended val- 
ues from various experiments and theories [61]. SE refers to 
semi-empirical calculations [53, 59, 60]. EXP and SE values 
are averaged over different j. 
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FIG. 2. (Color online) Pathways of 100 exemplary trajecto- 
ries of ionization dynamics of Xe at 4500 eV, 80 fs FWHM, 
and 5xl0 12 photons//im 2 . Multiphoton multiple ionization 
is described by a sequence of one-photon ionization (blue), 
Auger decay (green), and fluorescence (yellow). The gray 
background shows the Gaussian pulse profile of 80 fs FWHM. 



jectories of Xe at 4500 eV. The pulse duration is 80 fs 
full-width-at-half-maximum (FWHM) and the fluence is 
5xl0 12 photons//im 2 . The blue and green bars repre- 
sent photoionization and Auger (Coster-Kronig) decay, 
respectively, and the yellow dots indicate fluorescence. 
The bar colors are transparent, so darker colors mean 
that it is more probable to pass through those pathways. 
The ionization dynamics are obviously initiated by M- 
shell one-photon ionization, as shown by the blue bars 
between charge states zero and one. This is followed by 
a series of Auger decays, as shown by the green area 
above the initial blue bars, i.e., an Auger cascade after 
one-photon absorption [62]. Decay pathways after one- 
photon 3d-shell ionization of Xe have been studied exper- 
imentally [63, 64] and theoretically [65-67]. Note that the 
time scale of the Auger cascade ranges from 10 to 100 fs, 
as depicted in Fig. 2. Thus, more than one photon can 
be absorbed before the Auger cascade ends, as shown by 
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FIG. 3. (Color online) Charge state distribution of Xe at 
4500 eV as a function of the fluence. 
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FIG. 4. (Color online) Electron spectra of Xe at 4500 eV as a 
function of the fluence. The photoelectrons are above 1250 eV 
and the Auger (Coster-Kronig) electrons are below 1250 eV. 



sparse blue bars inside the green area, thus opening up 
new channels for ionization. As the charge state goes 
up, Auger decays become less likely [12], so photoioniza- 
tion becomes dominant for further ionization around the 
peak of the pulse profile. Fluorescence typically occurs 
at high charge states, when its rate overcomes the Auger- 
rate. The final charge states are formed in the middle of 
the latter half of the pulse. At the center of the pulse 
the charge states around +20 to +30 are formed, and at 
the end of the pulse they are distributed around +30 to 
+40. The pulse-weighted time-averaged charge state is 
+24 for this fluence case. 



C. Ion, electron, and photon spectra 

After calculating all 40 million atomic data parameters 
for one million configurations, rate equations for given 
XFEL parameters are solved in the Monte Carlo fash- 
ion described in Sec. II B. From the Monte Carlo simu- 
lations, we investigate ion, electron, and photon spectra 
by counting the final charge states, and photo-/Auger 
electrons and emitted photons in the energy bins. The 
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FIG. 5. (Color online) Fluorescence spectra of Xe at 4500 eV 
as a function of the fluence. The peak assignments are ex- 
plained in the text. 



pulse envelope is Gaussian and the pulse duration is 80 fs 
FWHM. In the regime of sequential ionization dynamics, 
the spectra are largely insensitive to the temporal pulse 
shape and the pulse duration [19]. The maximum fluence 
used is 10 13 photons//xm 2 . Even at the peak intensity of 
this fluence, the inverses of all photoionization rates are 
longer than 2.6 femtoseconds. The typical bandwidth of 
current XFEL sources operating in the hard x-ray regime 
is about 1% of the photon energy at LCLS [12, 13, 17] or 
tens of eV at SACLA in Japan [69]. If the XFEL band- 
width is given by 45 eV FWHM for a photon energy of 
4500 eV, then the coherence time is about 40 attoseconds. 
Thus the coherence time is much shorter than the time 
scale of the fastest photoionization process, which war- 
rants the use of the rate equation model in this regime. 

Figure 3 shows the charge state distribution of Xe at 
4500 eV after the x-ray pulse is over and all decay pro- 
cesses are completed. The vertical axis is the fluence 
varying from zero to 10 13 photons/zxm 2 , and the color 
indicates the fractional yield of a given charge state. The 
fractional yields are normalized such that the sum of all 
populations of ions and neutral atom is one. Therefore, 
one can observe that the total yield of ions is increasing 
as the fluence increases. For clarity, the population of 
neutral Xe is not shown in the plot. Near zero fluence, 
which corresponds to the synchrotron radiation limit, the 
charge state distribution peaks around +6 and +7. This 
distribution is due to the decay cascades of M-shell sin- 
gle vacancies [37, 68, 70]. When the fluence increases, 
the charge state distribution is shifted to higher charge 
states. The maximum charge state is +44, where all M-, 
N-, and O-shell electrons are ionized. When the pho- 
ton energy is not enough to ionize electrons by absorbing 
one photon, two-photon ionization may occur [13, 71] 
or resonantly excited states can play a role in generat- 
ing higher charge states [17]. As shown in Fig. 1, the 
4500-eV photon energy is large enough to ionize all elec- 
trons above the L-shell via one-photon absorption for all 
charge states. Therefore contributions from direct two- 
photon ionization and resonant pathways are negligible 



at 4500 eV. The maximum charge state of +44 can be 
reached via a sequence of one-photon processes. 



In Fig. 4, we plot electron spectra of Xe at 4500 eV, in- 
cluding both photoelectrons and Auger (Coster-Kronig) 
electrons. The vertical axis is the fluence and the hor- 
izontal axis is the electron kinetic energy spaced by 
the 10-eV width of the energy bins. Below ^0.5 keV, 
Auger (Coster-Kronig) electrons appear, and photoelec- 
trons have higher kinetic energy in the range from 1.5 
to 4.0 keV. There are three peak lines from photoelec- 
trons: 3.5, 3.6, and 3.8 keV, corresponding to photoion- 
ization from 3s, 3p, and 3c? of neutral Xe, respectively. 
From x-ray atomic data for neutral Xe [72], those lines 
are located at 3.4, 3.5, and 3.8 keV. Shake-off satellite 
structures in the electron spectra are not included in our 
calculations. The fine satellite structures shown in Fig. 4, 
especially in the photoelectron spectra, are due to differ- 
ent charge states and individual configurations. With 
increasing the charge state, the ionization threshold be- 
comes larger and accordingly the photoelectron kinetic 
energy becomes smaller. Therefore, as the fluence in- 
creases, higher charge states are formed and the photo- 
electron spectra extend to lower energies. 



Figure 5 shows fluorescence spectra of Xe at 4500 eV. 
Note that the width of the energy bins is 10 eV and the 
color bar is in the logarithmic scale. The comb structure 
comes from different charge states and the fringes are 
due to different electronic configurations. The photon 
spectra can be grouped according to different transition 
channels: i) strong lines below ~300 eV, ii) a plume from 
keV to 0.5 keV, iii) a comb structure plus a cloud of 
lines from 0.7 keV to 1.5 keV, and iv) comb lines from 
1.3 keV to 2.2 keV. To assign those parts, we plot in Fig. 6 
fluorescence energies for several transition channels as a 
function of the charge state. These energies are calcu- 
lated from orbital binding energies of the ground config- 
urations for given charge states. When the charge state 
is increased, transition energies between different shells 
(different quantum number n) are increased because the 
energy levels for n are approximately proportional to the 
square of the charge state. On the other hand, transi- 
tion energies between subshells of the same n (but differ- 
ent quantum number I) are decreased because electronic 
screening becomes less for higher charge states. There- 
fore, the four different groups in Fig. 5 can be assigned as 
follows: i) n=3 to n=3, ii) n—b to n=4, iii) n=4 to n=3, 
and iv) n=5 to n=3. Contrary to the electron spectra, 
the photon spectra extend to increasing energies as the 
fluence increases. If the energy resolution of the photon 
spectra is better than 10 eV as used in Fig. 5, then it 
is possible to observe these comb structures in the flu- 
orescence lines and to assign them to individual charge 
states, which was recently demonstrated at LCLS in the 
case of solid aluminum [18]. 



7 





1 /:* 










* * 
* » 
f > 








- 


\ ^ 








- 




w 












» » 
» ' 
» ' 












• » 
» » 
» • 












' M 






3p,d — 


3s,p 




II 




ob.p.u * is,p,u ------ 

i 4s.p,d -* 3s,p,d 




: 


'• : '• • i: 


: 5s.p,d -* 3s,p,d 

! i i 



0.5 1.0 1.5 2.0 2.5 

Photon energy (keV) 

FIG. 6. (Color online) Selected fluorescence energies of Xe as 
a function of the charge state. 

IV. CONCLUSION 

In this paper, we have implemented a Monte Carlo 
strategy for solving a rate equation model describing 
ionization dynamics induced by intense x-ray pulses. 
Photoionization cross sections, Auger (Coster-Kronig) 
rates, and fluorescence rates are calculated for all pos- 
sible multiple-hole configurations. Based on the pre- 
calculated table of all atomic data, Monte Carlo sam- 
pling finds probable pathways to reach the final charge 
states. Using the XATOM toolkit extended by this Monte 
Carlo method, we have investigated ionization dynam- 
ics of Xe in 4500-eV XFEL pulses. Detailed ionization 
and relaxation pathways have been depicted as a function 
of time. We have plotted the charge state distribution, 
photo-/Auger electron spectra, and fluorescence spectra 
as a function of fluence, whose range is experimentally ac- 
cessible. Near the upper end of this range, Xe at 4500 eV 
can be ionized up to +44 via a sequence of one-photon 



ionization and relaxation processes. 

Finally, we would like to briefly describe perspectives 
for further development. First, the current Monte Carlo 
implementation could be called "brute-force" since it cal- 
culates atomic data for all possible configurations and 
physical processes. The Monte Carlo sampling is applied 
for solving rate equations, but not for calculating atomic 
data. It is plausible to integrate both atomic data calcu- 
lation and search for probable pathways into the Monte 
Carlo procedure. In this way, atomic data are computed 
only when they are required. Second, it is important to 
include bound-to-bound photoexcitation processes in the 
model, which may play a crucial role in ionization dynam- 
ics at certain conditions. Because ionization thresholds 
have a broad range according to charge states and be- 
cause XFEL pulses typically have a broad bandwidth, 
resonance conditions may be easily satisfied. However, 
treatment of singly or multiply excited states is theoret- 
ically challenging, and moreover, inclusion of all these 
additional excited-state configurations is numerically de- 
manding. Third, relativistic effects are not considered in 
the present work, because M-, N-, and O-shell ioniza- 
tion of Xe may be described reasonably well within the 
HFS model. It is conceivable that spin-orbit energy split- 
tings open additional decay channels, but the configura- 
tional space will expand substantially when all splittings 
are taken into account. Work towards overcoming these 
challenges is in progress. 
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